This markdown document is meant to illustrate all the results presented in Ilboudo., et al (2018), all the preprocessing done on the data and the generation of tables as well as figures. The document is interactive and can be be easily navigated moving from one section to another.
Background: Sickle cell disease (SCD) is a monogenic disease caused by mutations in the β-globin gene. The complications related to the disease are systemic as they impact multiple organ systems. Our goal in this study was to identify metabolome changes contributing to SCD-related severity.
Employing both targeted and untargeted approaches, we profiled the plasma of 706 SCD patients using liquid chromatography tandem mass spectrometry. The cohort included 406 French patients of recent African descent and 300 African Americans from the southeastern US. In total, we measured the levels of 61 known and 2,100 unknown metabolites. We applied weighted gene correlation network analysis (WGCNA) algorithms to account for correlations among metabolites and identify specific metabolomic modules associated with SCD-related complications
Several chuncks of code for this analysis are adapted from WGCNA tutorial.
Load or install libraries necessary to perform analysis
source("http://bioconductor.org/biocLite.R")
biocLite(c("AnnotationDbi", "impute", "GO.db", "preprocessCore","WGCNA"))
install.packages(pkgs=c("WGCNA", "flashClust", "ggplot2",
"limma","factoextra","devtools","Biobase","sva","dplyr","magrittr"))
devtools::install_github("yilboudo/myUtils")
rm(list=ls())
options(stringsAsFactors = F)
library(myUtils)
library(limma)
library(WGCNA)
library(factoextra)
library(devtools)
library(Biobase)
library(sva)
library(dplyr)
library(gdata)
library(data.table)
allsamples_new <- read.table("GENMOD_OMG_metabolites_no_batch_correction.txt",sep="\t",h=T)
batches <- read.table("GENMOD_OMG_batch_information.txt",sep="\t",h=T)
#Load dplyr
slice <- dplyr::slice
#align data in the right order
x <- rownames(allsamples_new)
batches <- batches %>% slice(match(x, CLIENT_IDS))
batches <- as.data.frame(batches)
#image_name1 = 'all_cohorts_batches_pre_correction.pdf'
#pdf(file = image_name1, width = 10, height = 10)
res.pca_uncorrected <- prcomp(allsamples_new, scale = T)
fviz_pca_ind(res.pca_uncorrected, label="none", habillage=batches$Batches,
addEllipses=F, ellipse.level=0.95, palette = "igv") #+ xlim(-50,70) + ylim(-50,50)
#dev.off()
allsamples_new_combat <- run_combat(allsamples_new, batches, 2)
res.pca_corrected <- prcomp(allsamples_new_combat, scale = TRUE)
#save resulting output
allsamples_new_combat$CLIENT_IDs <- batches$CLIENT_IDS
write.table(allsamples_new_combat, file = "GENMOD_OMG_metabolites_with_batch_correction.txt",
quote = FALSE, row.names = FALSE, col.names = TRUE,sep="\t")
image_name2 = 'all_cohorts_batches_post_correction.pdf'
pdf(file = image_name2, width = 10, height = 10)
fviz_pca_ind(res.pca_corrected, label = "none", habillage=batches$Batches,
addEllipses=F, ellipse.level=0.95, palette = "igv") + xlim(-65,40) + ylim(-65,65)
dev.off()
#Load dataset
#Metabolites
allsamples_new_combat <- read.table("GENMOD_OMG_metabolites_with_batch_correction.txt",sep="\t",h=T,stringsAsFactors = F)
#Phenotypes
large_data_genmod_duke_n705 <-read.table("genmod_duke_phenotypes_n705_final.txt",sep="\t",h=T,stringsAsFactors = F)
#Clean metabolite data
rownames(allsamples_new_combat) <- allsamples_new_combat$CLIENT_IDs
allsamples_new_combat <- allsamples_new_combat[,-c(which(names(allsamples_new_combat)=="CLIENT_IDs"))]
#Apply Cyclic Loess Normalization
allsamples_new_combat_norm <- normalizeCyclicLoess(allsamples_new_combat, weights = NULL, span=0.7, iterations = 3, method = "fast")
allsamples_new_combat_norm_cl <- as.data.frame(allsamples_new_combat_norm)
#Apply clustering
sampleTree = hclust(dist(allsamples_new_combat_norm_cl), method = "average")
sizeGrWindow(12,9)
#pdf(file = "allsamples_Clustering.pdf", width = 12, height = 9)
#par(cex = 0.6);
#par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
cex.axis = 1.5, cex.main = 2)
abline(h = 52, col = "red")
#dev.off()
Re-cluster samples and remove outliers from dataset
# Determine cluster under the line
clust = cutreeStatic(sampleTree, cutHeight = 52, minSize = 10)
keepSamples = (clust==1)
allsamples_norm_cl_trimmed = allsamples_new_combat_norm_cl[keepSamples, ]
nMetabo = ncol(allsamples_norm_cl_trimmed)
nSamples = nrow(allsamples_norm_cl_trimmed)
Remove outliers from phenotype dataset
allTraits = large_data_genmod_duke_n705;
datTraits <- allTraits[match(rownames(allsamples_norm_cl_trimmed), rownames(allTraits)),]
#dim(datTraits)
collectGarbage();
Cluster samples showing phenotype information on blood traits and complications
sampleTree2 = hclust(dist(allsamples_norm_cl_trimmed), method = "average")
#Convert traits to a color representation: white means low, red means high, grey means missing entry
blood_traitColors = numbers2colors(datTraits[,c(which(names(datTraits)=="eGFR"),which(names(datTraits)=="Hematocrit"),
which(names(datTraits)=="HemoglobinF"),which(names(datTraits)=="MCH"),
which(names(datTraits)=="MCV"),which(names(datTraits)=="Platelets"),
which(names(datTraits)=="RBC"),which(names(datTraits)=="Reticulocyte"),
which(names(datTraits)=="Whitebloodcount"),which(names(datTraits)=="Bilirubin"),
which(names(datTraits)=="LDH"))],signed=T)
#blood_traitColors = suppressWarnings(numbers2colors(datTraits[,c(69,27,24)],signed=T))
# Convert traits to a color representation: white means low, red means high, grey means missing entry
complications_traitColors = numbers2colors(datTraits[,c(which(names(datTraits)=="Legulcer"),which(names(datTraits)=="Cholecystectomy"), which(names(datTraits)=="AscepticNecrosis"),which(names(datTraits)=="Priapism"),
which(names(datTraits)=="Retinopathy"),which(names(datTraits)=="Survival"))],signed=T)
Display graphics for blood traits and for complications
#pdf(file = "allsamples_with_bloodtraits_clustering2.pdf", width = 20, height = 15);
plotDendroAndColors(sampleTree2, blood_traitColors,
groupLabels = names(datTraits[,c(which(names(datTraits)=="eGFR"),which(names(datTraits)=="Hematocrit"),
which(names(datTraits)=="HemoglobinF"),which(names(datTraits)=="MCH"),
which(names(datTraits)=="MCV"),which(names(datTraits)=="Platelets"),
which(names(datTraits)=="RBC"),which(names(datTraits)=="Reticulocyte"),
which(names(datTraits)=="Whitebloodcount"),which(names(datTraits)=="Bilirubin"),
which(names(datTraits)=="LDH"))]), main = "Sample dendrogram and trait heatmap")
#dev.off()
#pdf(file = "allsamples_with_complication_traits_clustering.pdf", width = 20, height = 15);
plotDendroAndColors(sampleTree2, complications_traitColors,
groupLabels = names(datTraits[,c(which(names(datTraits)=="Legulcer"),
which(names(datTraits)=="Cholecystectomy"),
which(names(datTraits)=="AscepticNecrosis"),which(names(datTraits)=="Priapism"),
which(names(datTraits)=="Retinopathy"),which(names(datTraits)=="Survival"))]),
main = "Sample dendrogram and trait heatmap")
#dev.off()
Save resulting data frames
save(allsamples_norm_cl_trimmed, datTraits, file = "allsamples-01-dataInput.RData")
powers = c(c(1:10), seq(from = 12, to=20, by=2))
sft = pickSoftThreshold(allsamples_norm_cl_trimmed,corFnc = cor, networkType = 'signed',powerVector = powers, verbose = 5)
Choosing the soft-thresholding power: analysis of network topology
#pdf(file = "allsamples_power.pdf", width = 12, height = 9)
par(mfrow = c(1,2));
cex1 = 0.8;
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=cex1,col="red");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.80,col="red")
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
#dev.off()
#Co-expression similarity and adjacency
softPower = 6;
adjacency = adjacency(allsamples_norm_cl_trimmed, power = softPower);
#Topological Overlap Matrix (TOM)
TOM = TOMsimilarity(adjacency);
dissTOM = 1-TOM
#Clustering using TOM
geneTree = hclust(as.dist(dissTOM), method = "average");
# Plot the resulting clustering tree (dendrogram)
pdf(file = "allsamples_metabo_networks.pdf", width = 12, height = 9)
plot(geneTree, xlab="", sub="", main = "Metabolite clustering on TOM-based dissimilarity",
labels = FALSE, hang = 0.04);
dev.off()
Module identification using dynamic tree cut
minModuleSize = 7;
# Module identification using dynamic tree cut:
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
deepSplit = 2, pamRespectsDendro = FALSE,
minClusterSize = minModuleSize);
table(dynamicMods)
# Convert numeric lables into colors
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)
Plot Dendogram with colors underneath
# Plot the dendrogram and colors underneath
sizeGrWindow(8,6)
pdf(file = "allsamples_metabo_networks_color.pdf", width = 12, height = 9)
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05,
main = "Metabolite dendrogram and module colors")
dev.off()
lnames3 = load(file = "~/Desktop/Metabolite_All_Data_dec52017/module_preservation_stats.RData")
preserved_modules <- preservations_stats_all_df[preservations_stats_all_df$V2 > 10,]
restGenes= (moduleColors %in% preserved_modules$V1) & !(moduleColors %in% "grey") & !(moduleColors %in% "gold")
SubGeneNames <- colnames(allsamples_norm_cl_trimmed)
dissTOM = 1-TOMsimilarityFromExpr(allsamples_norm_cl_trimmed[,restGenes], power = 6);
# Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
plotTOM = dissTOM^10;
colnames(dissTOM) =rownames(dissTOM) =SubGeneNames[restGenes]
hier1=hclust(as.dist(dissTOM), method="average" )
# Set diagonal to NA for a nicer plot
diag(plotTOM) = NA;
# Call the plot function
#sizeGrWindow(9,9)
#pdf(file = "all_samples_network_plot2.pdf", width = 10, height = 10)
TOMplot(plotTOM, hier1, moduleColors[restGenes], main = "Network heatmap plot, all genes")
#dev.off()
MEList = moduleEigengenes(allsamples_norm_cl_trimmed[,restGenes], colors = moduleColors[restGenes])
MEs = MEList$eigengenes
MET=orderMEs(MEs)
#pdf(file = "all_samples_eigengenes_plot2.pdf", width = 10, height = 10)
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2))
#dev.off()
lnames1 = load(file = "allsamples-01-dataInput.RData");
#The variable lnames contains the names of loaded variables.
lnames1
# Load network data saved in the second part.
lnames2 = load(file = "allsamples-02-networkConstruction-stepByStep.RData");
datTraits_for_assoc <- datTraits
# Define numbers of metabolites and samples
nGenes = ncol(allsamples_norm_cl_trimmed);
nSamples = nrow(allsamples_norm_cl_trimmed);
# Recalculate MEs with color labels
MEs0 = moduleEigengenes(allsamples_norm_cl_trimmed, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
rownames(MEs) <- rownames(allsamples_norm_cl_trimmed)
#inverse normal transform eigengenes
MEs_for_assoc <- as.data.frame(apply(MEs,2,invnorm))
#Define variable for blood traits, and covariates
blood_traits1 = "eGFR"
cov_names1="SCD_geno,HU"
blood_traits = c("Hematocrit","HemoglobinF","MCH","MCV","MCHC","Platelets","RBC","Reticulocyte","Whitebloodcount", "Bilirubin", "LDH")
cov_names="SCD_geno,HU,Age,Sex"
#Create matrices to store summary statistics values
blood_traits_associations_pval <- matrix(NA,dim(MEs)[2],length(blood_traits))
blood_traits_associations_effect <- matrix(NA,dim(MEs)[2],length(blood_traits))
blood_traits_associations_n <- matrix(NA,dim(MEs)[2],length(blood_traits))
blood_traits_associations_pval1 <- matrix(NA,dim(MEs)[2],length(blood_traits1))
blood_traits_associations_effect1 <- matrix(NA,dim(MEs)[2],length(blood_traits1))
blood_traits_associations_n1 <- matrix(NA,dim(MEs)[2],length(blood_traits1))
#Run association for most traits except eGFR because covariates are not the same
for (i in 1:length(blood_traits)) {
datTraits_trait <- datTraits_for_assoc[complete.cases(datTraits_for_assoc[,which(colnames(datTraits_for_assoc)==blood_traits[i])]),]
keeping_indiv <- rownames(datTraits_trait)
datTraits_trait$Zpheno <- invnorm(datTraits_trait[,which(names(datTraits_trait)==blood_traits[i])])
MEsTrait <- MEs_for_assoc[(rownames(MEs_for_assoc) %in% keeping_indiv),]
covariates <- strsplit(cov_names, split=",")[[1]]
cov=datTraits_trait[,covariates]
observed_pval <- sapply(MEsTrait, linearRegcov_p, y=datTraits_trait$Zpheno, data=cov)
std_err <- sapply(MEsTrait, linearRegcov_stderr, y=datTraits_trait$Zpheno, data=cov)
effect_size <- sapply(MEsTrait, linearRegcov_effect, y=datTraits_trait$Zpheno, data=cov)
sample_size <- sapply(MEsTrait, linearRegcov_samplesize, y=datTraits_trait$Zpheno, data=cov)
blood_traits_associations_pval[,i] <- observed_pval
blood_traits_associations_effect[,i] <- effect_size
blood_traits_associations_n[,i] <- sample_size
}
#Run association for eGFR
for (i in 1:length(blood_traits1)) {
datTraits_trait <- datTraits_for_assoc[complete.cases(datTraits_for_assoc[,which(colnames(datTraits_for_assoc)==blood_traits1[i])]),]
keeping_indiv <- rownames(datTraits_trait)
datTraits_trait$Zpheno <- invnorm(datTraits_trait[,which(names(datTraits_trait)==blood_traits1[i])])
MEsTrait <- MEs_for_assoc[(rownames(MEs_for_assoc) %in% keeping_indiv),]
covariates <- strsplit(cov_names1, split=",")[[1]]
cov=datTraits_trait[,covariates]
observed_pval <- sapply(MEsTrait, linearRegcov_p, y=datTraits_trait$Zpheno, data=cov)
std_err <- sapply(MEsTrait, linearRegcov_stderr, y=datTraits_trait$Zpheno, data=cov)
effect_size <- sapply(MEsTrait, linearRegcov_effect, y=datTraits_trait$Zpheno, data=cov)
sample_size <- sapply(MEsTrait, linearRegcov_samplesize, y=datTraits_trait$Zpheno, data=cov)
blood_traits_associations_pval1[,i] <- observed_pval
blood_traits_associations_effect1[,i] <- effect_size
blood_traits_associations_n1[,i] <- sample_size
}
blood_traits_associations_pval <- cbind(blood_traits_associations_pval, blood_traits_associations_pval1)
blood_traits_associations_effect <- cbind(blood_traits_associations_effect, blood_traits_associations_effect1)
blood_traits_associations_n <- cbind(blood_traits_associations_n, blood_traits_associations_n1)
blood_traits_associations <- as.data.frame(blood_traits_associations_pval)
names(blood_traits_associations) <- c(blood_traits , "eGFR")
rownames(blood_traits_associations) <- names(MEs)
blood_traits_associations_cor <- as.data.frame(blood_traits_associations_effect)
names(blood_traits_associations_cor) <- c(blood_traits , "eGFR")
rownames(blood_traits_associations_cor) <- names(MEs)
#Associations with complications
complications = c("Legulcer","Cholecystectomy","AscepticNecrosis","Priapism", "Retinopathy","Survival")
complications_associations_pval <- matrix(NA,dim(MEs)[2],length(complications))
complications_associations_effect <- matrix(NA,dim(MEs)[2],length(complications))
complications_associations_n <- matrix(NA,dim(MEs)[2],length(complications))
for (i in 1:length(complications)) {
datTraits_trait <- datTraits_for_assoc[complete.cases(datTraits_for_assoc[,which(colnames(datTraits_for_assoc)==complications[i])]),]
keeping_indiv <- rownames(datTraits_trait)
datTraits_trait$Zpheno <- datTraits_trait[,which(names(datTraits_trait)==complications[i])]
MEsTrait <- MEs_for_assoc[(rownames(MEs_for_assoc) %in% keeping_indiv),]
covariates <- strsplit(cov_names, split=",")[[1]]
cov=datTraits_trait[,covariates]
observed_pval <- sapply(MEsTrait, logisticRegcov_p, y=datTraits_trait$Zpheno, data=cov)
effect_size <- sapply(MEsTrait, logisticRegcov_effect, y=datTraits_trait$Zpheno, data=cov)
sample_size <- sapply(MEsTrait, logisticRegcov_samplesize, y=datTraits_trait$Zpheno, data=cov)
complications_associations_pval[,i] <- observed_pval
complications_associations_effect[,i] <- effect_size
complications_associations_n[,i] <- sample_size
}
complications_associations <- as.data.frame(complications_associations_pval)
names(complications_associations) <- complications
rownames(complications_associations) <- names(MEs)
complications_associations_cor <- as.data.frame(complications_associations_effect)
names(complications_associations_cor) <- complications
rownames(complications_associations_cor) <- names(MEs)
pvals <- cbind(complications_associations,blood_traits_associations)
cor <- cbind(complications_associations_cor,blood_traits_associations_cor)
Display associations statistics of modules and blood traits and complications in a heatmap
png(file="allsamples_metabo_traits_assoc_blood_final.png", units = "px",width=9000, height=9000, res=350, pointsize=15)
# Will display correlations and their p-values
textMatrix = paste(signif(as.matrix(pvals), 2), "\n(",
signif(as.matrix(cor), 1), ")", sep = "");
dim(textMatrix) = dim(cor)
par(mar = c(7, 8.5, 2, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = cor[,which(names(pvals)=="Hematocrit"):which(names(pvals)=="eGFR")],
xLabels = names(pvals)[which(names(pvals)=="Hematocrit"):which(names(pvals)=="eGFR")],
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = F,
colors = blueWhiteRed(50),
textMatrix = textMatrix[,which(names(pvals)=="Hematocrit"):which(names(pvals)=="eGFR")],
setStdMargins = F,
cex.text = 0.7,
zlim = c(-0.60, 0.40),
main = paste("Module-trait relationships"))
dev.off()
png(file="allsamples_metabo_traits_assoc_complications_final.png", units = "px",width=9000, height=9000, res=350, pointsize=15)
# Will display correlations and their p-values
textMatrix = paste(signif(as.matrix(pvals), 2), "\n(",
signif(as.matrix(cor), 1), ")", sep = "");
dim(textMatrix) = dim(cor)
par(mar = c(7, 8.5, 2, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = cor[,which(names(pvals)=="Legulcer"):which(names(pvals)=="Survival")], #moduleTraitCor[,-c(1,15:19)],
xLabels = names(pvals)[which(names(pvals)=="Legulcer"):which(names(pvals)=="Survival")],
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = F,
colors = blueWhiteRed(50),
textMatrix = textMatrix[,which(names(pvals)=="Legulcer"):which(names(pvals)=="Survival")],
setStdMargins = F,
cex.text = 0.7,
zlim = c(0,2),
main = paste("Module-trait relationships"))
dev.off()
Find top module per trait
#Top module with strong traits associations
topmodule <- rownames(pvals[pvals$Legulcer< (0.05/(66*15)),])
topmodule <- append(topmodule,rownames(pvals[pvals$Cholecystectomy< (0.05/(66*15)),]))
topmodule <- append(topmodule,rownames(pvals[pvals$AscepticNecrosis< (0.05/(66*15)),]))
topmodule <- append(topmodule,rownames(pvals[pvals$Priapism< (0.05/(66*15)),]))
topmodule <- append(topmodule,rownames(pvals[pvals$Retinopathy< (0.05/(66*15)),]))
topmodule <- append(topmodule,rownames(pvals[pvals$Hematocrit< (0.05/(66*15)),]))
topmodule <- append(topmodule,rownames(pvals[pvals$HemoglobinF< (0.05/(66*15)),]))
topmodule <- append(topmodule,rownames(pvals[pvals$MCH< (0.05/(66*15)),]))
topmodule <- append(topmodule,rownames(pvals[pvals$MCV< (0.05/(66*15)),]))
topmodule <- append(topmodule,rownames(pvals[pvals$Platelets< (0.05/(66*15)),]))
topmodule <- append(topmodule,rownames(pvals[pvals$RBC< (0.05/(66*15)),]))
topmodule <- append(topmodule,rownames(pvals[pvals$Reticulocyte< (0.05/(66*15)),]))
topmodule <- append(topmodule,rownames(pvals[pvals$Whitebloodcount< (0.05/(66*15)),]))
topmodule <- append(topmodule,rownames(pvals[pvals$Survival< (0.05/(66*15)),]))
topmodule <- append(topmodule,rownames(pvals[pvals$eGFR< (0.05/(66*15)),]))
table(topmodule)
modNames = substring(names(MEs), 3)
eigens <- names(MEs)
dim(datTraits)
#Metabolites to module eigenvalues
metabo_modulemembership_cor = matrix(NA,length(1:length(allsamples_norm_cl_trimmed)),dim(MEs_for_assoc)[2])
metabo_modulemembership_pval = matrix(NA,length(1:length(allsamples_norm_cl_trimmed)),dim(MEs_for_assoc)[2])
#
datTraits_trait <- datTraits_for_assoc[complete.cases(datTraits_for_assoc[,which(colnames(datTraits_for_assoc)=="Sex")]),]
#
covariates <- strsplit(cov_names, split=",")[[1]]
cov=datTraits_trait[,covariates]
#
for (i in 1:length(eigens)) {
observed_pval <- sapply(allsamples_norm_cl_trimmed, linearRegcov_p, y=MEs_for_assoc[,eigens[i]], data=cov)
effect_size <- sapply(allsamples_norm_cl_trimmed, linearRegcov_effect, y=MEs_for_assoc[,eigens[i]], data=cov)
metabo_modulemembership_pval[,i] <- observed_pval
metabo_modulemembership_cor[,i] <- effect_size
}
metabo_modulemembership_pval_final <- as.data.frame(metabo_modulemembership_pval)
names(metabo_modulemembership_pval_final) <- paste(names(MEs_for_assoc),".pval",sep="")
rownames(metabo_modulemembership_pval_final) <- colnames(allsamples_norm_cl_trimmed)
metabo_modulemembership_cor_final <- as.data.frame(metabo_modulemembership_cor)
names(metabo_modulemembership_cor_final) <- paste(names(MEs_for_assoc),".beta",sep="")
rownames(metabo_modulemembership_cor_final) <- colnames(allsamples_norm_cl_trimmed)
geneModuleMembership=metabo_modulemembership_cor_final
MMPvalue = metabo_modulemembership_pval_final
save(eigens,modNames,metabo_modulemembership_cor, metabo_modulemembership_pval,
metabo_modulemembership_pval_final,metabo_modulemembership_cor_final,
file = "Module_Memberships_Associations.RData")
load("Module_Memberships_Associations.RData")
#Metabolites to traits
metabo_bloodtrait_cor = matrix(NA,length(1:length(allsamples_norm_cl_trimmed)), length(blood_traits))
metabo_bloodtrait_pval = matrix(NA,length(1:length(allsamples_norm_cl_trimmed)), length(blood_traits))
metabo_bloodtrait_n = matrix(NA,length(1:length(allsamples_norm_cl_trimmed)), length(blood_traits))
metabo_bloodtrait_cor1 = matrix(NA,length(1:length(allsamples_norm_cl_trimmed)), length(blood_traits1))
metabo_bloodtrait_pval1 = matrix(NA,length(1:length(allsamples_norm_cl_trimmed)), length(blood_traits1))
metabo_bloodtrait_n1 = matrix(NA,length(1:length(allsamples_norm_cl_trimmed)), length(blood_traits1))
for (i in 1:length(blood_traits)) {
datTraits_trait <- datTraits_for_assoc[complete.cases(datTraits_for_assoc[,which(colnames(datTraits_for_assoc)==blood_traits[i])]),]
keeping_indiv <- rownames(datTraits_trait)
datTraits_trait$Zpheno <- invnorm(datTraits_trait[,which(names(datTraits_trait)==blood_traits[i])])
allsamples_norm_cl_trimmedTrait <- allsamples_norm_cl_trimmed[(rownames(allsamples_norm_cl_trimmed) %in% keeping_indiv),]
covariates <- strsplit(cov_names, split=",")[[1]]
cov=datTraits_trait[,covariates]
observed_pval <- sapply(allsamples_norm_cl_trimmedTrait, linearRegcov_p, y=datTraits_trait$Zpheno, data=cov)
std_err <- sapply(allsamples_norm_cl_trimmedTrait, linearRegcov_stderr, y=datTraits_trait$Zpheno, data=cov)
effect_size <- sapply(allsamples_norm_cl_trimmedTrait, linearRegcov_effect, y=datTraits_trait$Zpheno, data=cov)
sample_size <- sapply(allsamples_norm_cl_trimmedTrait, linearRegcov_samplesize, y=datTraits_trait$Zpheno, data=cov)
metabo_bloodtrait_pval [,i] <- observed_pval
metabo_bloodtrait_cor[,i] <- effect_size
metabo_bloodtrait_n[,i] <- sample_size
}
for (i in 1:length(blood_traits1)) {
datTraits_trait <- datTraits_for_assoc[complete.cases(datTraits_for_assoc[,which(colnames(datTraits_for_assoc)==blood_traits1[i])]),]
keeping_indiv <- rownames(datTraits_trait)
datTraits_trait$Zpheno <- invnorm(datTraits_trait[,which(names(datTraits_trait)==blood_traits1[i])])
allsamples_norm_cl_trimmedTrait <- allsamples_norm_cl_trimmed[(rownames(allsamples_norm_cl_trimmed) %in% keeping_indiv),]
covariates <- strsplit(cov_names1, split=",")[[1]]
cov=datTraits_trait[,covariates]
observed_pval <- sapply(allsamples_norm_cl_trimmedTrait, linearRegcov_p, y=datTraits_trait$Zpheno, data=cov)
std_err <- sapply(allsamples_norm_cl_trimmedTrait, linearRegcov_stderr, y=datTraits_trait$Zpheno, data=cov)
effect_size <- sapply(allsamples_norm_cl_trimmedTrait, linearRegcov_effect, y=datTraits_trait$Zpheno, data=cov)
sample_size <- sapply(allsamples_norm_cl_trimmedTrait, linearRegcov_samplesize, y=datTraits_trait$Zpheno, data=cov)
metabo_bloodtrait_pval1 [,i] <- observed_pval
metabo_bloodtrait_cor1[,i] <- effect_size
metabo_bloodtrait_n1[,i] <- sample_size
}
metabo_bloodtrait_cor <- cbind(metabo_bloodtrait_cor,metabo_bloodtrait_cor1)
metabo_bloodtrait_cor_final <- as.data.frame(metabo_bloodtrait_cor)
names(metabo_bloodtrait_cor_final) <- paste(c(blood_traits,"eGFR"),".beta",sep="")
rownames(metabo_bloodtrait_cor_final) <- colnames(allsamples_norm_cl_trimmed)
metabo_bloodtrait_pval <- cbind(metabo_bloodtrait_pval,metabo_bloodtrait_pval1)
metabo_bloodtrait_pval_final <- as.data.frame(metabo_bloodtrait_pval)
names(metabo_bloodtrait_pval_final) <- paste(c(blood_traits,"eGFR"),".pval",sep="")
rownames(metabo_bloodtrait_pval_final) <- colnames(allsamples_norm_cl_trimmed)
metabo_bloodtrait_n <- cbind(metabo_bloodtrait_n,metabo_bloodtrait_n1)
metabo_bloodtrait_n_final <- as.data.frame(metabo_bloodtrait_n)
names(metabo_bloodtrait_n_final) <- paste(c(blood_traits,"eGFR"),".N",sep="")
rownames(metabo_bloodtrait_n_final) <- colnames(allsamples_norm_cl_trimmed)
save(metabo_bloodtrait_n_final, metabo_bloodtrait_pval_final, metabo_bloodtrait_cor_final, file = "Module_Blood_Traits_Associations.RData")
#Metabolites to complications
metabo_complications_cor = matrix(NA,length(1:length(allsamples_norm_cl_trimmed)), length(complications))
metabo_complications_pval = matrix(NA,length(1:length(allsamples_norm_cl_trimmed)), length(complications))
metabo_complications_n = matrix(NA,length(1:length(allsamples_norm_cl_trimmed)), length(complications))
for (i in 1:length(complications)) {
datTraits_trait <- datTraits_for_assoc[complete.cases(datTraits_for_assoc[,which(colnames(datTraits_for_assoc)==complications[i])]),]
keeping_indiv <- rownames(datTraits_trait)
datTraits_trait$Zpheno <- datTraits_trait[,which(names(datTraits_trait)==complications[i])]
allsamples_norm_cl_trimmedTrait <- allsamples_norm_cl_trimmed[(rownames(allsamples_norm_cl_trimmed) %in% keeping_indiv),]
covariates <- strsplit(cov_names, split=",")[[1]]
cov=datTraits_trait[,covariates]
observed_pval <- sapply(allsamples_norm_cl_trimmedTrait, logisticRegcov_p, y=datTraits_trait$Zpheno, data=cov)
effect_size <- sapply(allsamples_norm_cl_trimmedTrait, logisticRegcov_effect, y=datTraits_trait$Zpheno, data=cov)
sample_size <- sapply(allsamples_norm_cl_trimmedTrait, logisticRegcov_samplesize, y=datTraits_trait$Zpheno, data=cov)
metabo_complications_pval [,i] <- observed_pval
metabo_complications_cor[,i] <- effect_size
metabo_complications_n[,i] <- sample_size
}
metabo_complications_cor_final <- as.data.frame(metabo_complications_cor)
names(metabo_complications_cor_final) <- paste(complications,".beta",sep="")
rownames(metabo_complications_cor_final) <- colnames(allsamples_norm_cl_trimmed)
metabo_complications_pval_final <- as.data.frame(metabo_complications_pval)
names(metabo_complications_pval_final) <- paste(complications,".pval",sep="")
rownames(metabo_complications_pval_final) <- colnames(allsamples_norm_cl_trimmed)
metabo_complications_n_final <- as.data.frame(metabo_complications_n)
names(metabo_complications_n_final) <- paste(complications,".N",sep="")
rownames(metabo_complications_n_final) <- colnames(allsamples_norm_cl_trimmed)
save(metabo_complications_n_final, metabo_complications_pval_final, metabo_complications_cor_final, file = "Module_Complications_Associations.RData")
lcomp <- load("Module_Complications_Associations.RData")
metabo_annotation <- fread("~/Desktop/Metabolite_All_Data_dec52017/wgcna_scd/data/full_list_from_analysis_annotated_FINAL.txt",sep="\t",h=F,stringsAsFactors = F)
names(metabo_annotation) <- c("metabo_name", "hmdb_id_input", "metabo_super_class","metabo_class","metabo_subclass",
"metabo_biofunction","metabo_accession_secondary","metabo_disease","metabo_gene",
"metabo_description")
metabo_bloodtrait_pval_final$metabo_name <- rownames(metabo_bloodtrait_pval_final)
metabo_complications_pval_final$metabo_name <- rownames(metabo_complications_pval_final)
metabo_modulemembership_pval_final$metabo_name <- rownames(metabo_modulemembership_pval_final)
metabo_bloodtrait_pval_final_annotation <- merge(metabo_bloodtrait_pval_final, metabo_annotation, by="metabo_name")
metabo_complications_pval_final_annotation <- merge(metabo_complications_pval_final, metabo_annotation, by="metabo_name")
metabo_modulemembership_pval_final_annotation <- merge(metabo_modulemembership_pval_final, metabo_annotation, by="metabo_name")
metabo_bloodtrait_cor_final$metabo_name <- rownames(metabo_bloodtrait_cor_final)
metabo_complications_cor_final$metabo_name <- rownames(metabo_complications_cor_final)
metabo_modulemembership_cor_final$metabo_name <- rownames(metabo_modulemembership_cor_final)
metabo_bloodtrait_cor_pval_final_annotation <- merge(metabo_bloodtrait_cor_final, metabo_bloodtrait_pval_final_annotation, by="metabo_name")
metabo_complications_cor_pval_final_annotation <- merge(metabo_complications_cor_final, metabo_complications_pval_final_annotation, by="metabo_name")
metabo_modulemembership_cor_pval_final_annotation <- merge(metabo_modulemembership_cor_final, metabo_modulemembership_pval_final_annotation, by="metabo_name")
metabolite_group_order <- colnames(allsamples_norm_cl_trimmed)
metabolite_color_match <- cbind(metabolite_group_order,moduleColors)
#########
metabo_modulemembership_cor_pval_final_annotation <- metabo_modulemembership_cor_pval_final_annotation[match(metabolite_color_match[,1], metabo_modulemembership_cor_pval_final_annotation$metabo_name),]
metabo_bloodtrait_cor_pval_final_annotation <- metabo_bloodtrait_cor_pval_final_annotation[match(metabolite_color_match[,1], metabo_bloodtrait_cor_pval_final_annotation$metabo_name),]
metabo_complications_cor_pval_final_annotation <- metabo_complications_cor_pval_final_annotation[match(metabolite_color_match[,1], metabo_complications_cor_pval_final_annotation$metabo_name),]
# Create the starting data frame
geneInfo0 = data.frame(
moduleColor = moduleColors,
metabo_bloodtrait_cor_pval_final_annotation)
geneInfo1 = data.frame(
moduleColor = moduleColors,
metabo_complications_cor_pval_final_annotation)
geneInfo2 = data.frame(
moduleColor = moduleColors,
metabo_modulemembership_cor_pval_final_annotation)
save(metabo_modulemembership_cor_pval_final_annotation, metabo_bloodtrait_cor_pval_final_annotation,
metabo_complications_cor_pval_final_annotation, geneInfo0,geneInfo1,geneInfo2,file = "allsamples-03-associations_annotation_external_traits.RData")
write.csv(geneInfo0, file = "metabo_bloodtrait_Info_FINAL2.csv",row.names = F)
write.csv(geneInfo1, file = "metabo_complication_Info_FINAL2.csv",row.names = F)
write.csv(geneInfo2, file = "metabo_module_Info_FINAL2.csv",row.names = F)
no.indivs = dim(allsamples_norm_cl_trimmed)[[1]]
##Plotting results
set.seed(1)
pdf(file = "Robustness_dropping_samples.pdf", width = 12, height = 9)
par(mfrow=c(8,2), mar=c(2,2,2,1))
for (i in 1:8) {
subsetindivs=sample(1:no.indivs,344,replace=F)
robustness_drop_samples(metabolites_data = allsamples_norm_cl_trimmed[subsetindivs,],
module_colors = moduleColors,power1=6 ,main1="SCD Metabolites Network",
filename=paste("metabolites_clust_",i,".txt",sep=""))
}
dev.off()
set.seed(1)
for (i in 1:1) {
subsetindivs=sample(1:no.indivs,344,replace=F)
multiExpr = list(original = list(data = allsamples_norm_cl_trimmed), dropped_sample = list(data = allsamples_norm_cl_trimmed[subsetindivs,]))
multiColor = list(original = moduleColors)
system.time( {
mp = modulePreservation(multiExpr, multiColor,
referenceNetworks = 1,
networkType = "signed",
nPermutations = 200,
randomSeed = 1,
quickCor = 0,
verbose = 2)
} );
save(mp, file = paste("modulePreservation.RData",sep="",i));
ref = 1
test = 2
statsObs = cbind(mp$quality$observed[[ref]][[test]][, -1], mp$preservation$observed[[ref]][[test]][, -1])
statsZ = cbind(mp$quality$Z[[ref]][[test]][, -1], mp$preservation$Z[[ref]][[test]][, -1]);
assign(paste("preservations_stats",i,sep=""), as.data.frame(cbind(statsObs[, c("medianRank.pres", "medianRank.qual")],
signif(statsZ[, c("Zsummary.pres", "Zsummary.qual")], 2))))
}
write.table(preservations_stats1 ,file= "preservations_stats1.txt",quote = FALSE, row.names = T, col.names = TRUE,sep="\t")
write.table(preservations_stats2 ,file= "preservations_stats2.txt",quote = FALSE, row.names = T, col.names = TRUE,sep="\t")
write.table(preservations_stats3 ,file= "preservations_stats3.txt",quote = FALSE, row.names = T, col.names = TRUE,sep="\t")
write.table(preservations_stats4 ,file= "preservations_stats4.txt",quote = FALSE, row.names = T, col.names = TRUE,sep="\t")
write.table(preservations_stats5 ,file= "preservations_stats5.txt",quote = FALSE, row.names = T, col.names = TRUE,sep="\t")
write.table(preservations_stats6 ,file= "preservations_stats6.txt",quote = FALSE, row.names = T, col.names = TRUE,sep="\t")
write.table(preservations_stats7 ,file= "preservations_stats7.txt",quote = FALSE, row.names = T, col.names = TRUE,sep="\t")
write.table(preservations_stats8 ,file= "preservations_stats8.txt",quote = FALSE, row.names = T, col.names = TRUE,sep="\t")
preservations_stats_all <- cbind(preservations_stats1$Zsummary.pres +preservations_stats2$Zsummary.pres +
preservations_stats3$Zsummary.pres +
preservations_stats4$Zsummary.pres +
preservations_stats5$Zsummary.pres +
preservations_stats6$Zsummary.pres +
preservations_stats7$Zsummary.pres +
preservations_stats8$Zsummary.pres)/8
preservations_stats_all_df <- as.data.frame(cbind(rownames(preservations_stats1),preservations_stats_all))
preservations_stats_all_df$V2 <- as.numeric(preservations_stats_all_df$V2)
preserved_modules <- preservations_stats_all_df[preservations_stats_all_df$V2 > 10,]
preserved_not_modules <- preservations_stats_all_df[preservations_stats_all_df$V2 < 10,]
save(preservations_stats_all_df, file = "module_preservation_stats.RData")
png(file="TopModule_IntramodularStrength_Complications.png", units = "px",width=3000, height=2000, res=350, pointsize=9)
scatterplots_metabo(module_color="lightcyan1",trait="Cholecystectomy.beta",metabolite_significance=metabo_complications_cor_final,
module_membership_cor=metabo_modulemembership_cor_final,alt_color="blue")
#scatterplots_metabo(module_color="coral1",trait="Legulcer.beta",dataset=metabo_complications_cor_final,alt_color="")
dev.off()
png(file="TopModule_IntramodularStrength_Blood_Hematocrit.png", units = "px",width=3000, height=2000, res=350, pointsize=9)
par(mfrow=c(3,5))
scatterplots_metabo(module_color="darkorange",trait="Hematocrit.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="yellow4",trait="Hematocrit.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="lightcyan",trait="Hematocrit.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="mediumorchid",trait="Hematocrit.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="firebrick4",trait="Hematocrit.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="maroon",trait="Hematocrit.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkgreen",trait="Hematocrit.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkslateblue",trait="Hematocrit.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="magenta",trait="Hematocrit.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="blue",trait="Hematocrit.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="brown4",trait="Hematocrit.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="skyblue1",trait="Hematocrit.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkgrey",trait="Hematocrit.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="lightcoral",trait="Hematocrit.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="sienna3",trait="Hematocrit.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
dev.off()
png(file="TopModule_IntramodularStrength_Blood_HemoglobinF.png", units = "px",width=3000, height=2000, res=350, pointsize=9)
scatterplots_metabo(module_color="darkslateblue",trait="HemoglobinF.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
dev.off()
png(file="TopModule_IntramodularStrength_Blood_RBC.png", units = "px",width=3000, height=2000, res=350, pointsize=9)
par(mfrow=c(3,4))
scatterplots_metabo(module_color="darkorange",trait="RBC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="yellow4",trait="RBC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="lightcyan",trait="RBC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="mediumorchid",trait="RBC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="firebrick4",trait="RBC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkolivegreen4",trait="RBC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="maroon",trait="RBC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkslateblue",trait="RBC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="skyblue1",trait="RBC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkgrey",trait="RBC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="lightcoral",trait="RBC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="sienna3",trait="RBC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
dev.off()
png(file="TopModule_IntramodularStrength_Blood_eGFR.png", units = "px",width=3000, height=2000, res=350, pointsize=9)
par(mfrow=c(3,7))
scatterplots_metabo(module_color="plum",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="greenyellow",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="honeydew1",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="skyblue3",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkorange",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="royalblue",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="lightcyan",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="blue")
scatterplots_metabo(module_color="darkolivegreen4",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="cyan",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="grey60",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="maroon",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkgreen",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkslateblue",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="magenta",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="coral1",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkmagenta",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="blue",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="coral2",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="brown",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="plum2",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="grey",trait="eGFR.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
dev.off()
png(file="TopModule_IntramodularStrength_Blood_MCH.png", units = "px",width=3000, height=2000, res=350, pointsize=9)
scatterplots_metabo(module_color="lightcyan",trait="MCH.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="blue")
dev.off()
png(file="TopModule_IntramodularStrength_Blood_MCHC.png", units = "px",width=3000, height=2000, res=350, pointsize=9)
par(mfrow=c(3,4))
scatterplots_metabo(module_color="yellow4",trait="MCHC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="lightcyan",trait="MCHC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="blue")
scatterplots_metabo(module_color="mediumorchid",trait="MCHC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="lightsteelblue1",trait="MCHC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="skyblue1",trait="MCHC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkgrey",trait="MCHC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="mediumpurple2",trait="MCHC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkorange2",trait="MCHC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="lightcoral",trait="MCHC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkred",trait="MCHC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="sienna3",trait="MCHC.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
dev.off()
png(file="TopModule_IntramodularStrength_Blood_MCV.png", units = "px",width=3000, height=2000, res=350, pointsize=9)
scatterplots_metabo(module_color="lightcyan",trait="MCV.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="blue")
dev.off()
png(file="TopModule_IntramodularStrength_Blood_Platelets.png", units = "px",width=3000, height=2000, res=350, pointsize=9)
par(mfrow=c(2,3))
scatterplots_metabo(module_color="darkorange2",trait="Platelets.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="skyblue",trait="Platelets.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkviolet",trait="Platelets.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
dev.off()
png(file="TopModule_IntramodularStrength_Blood_Reticulocyte.png", units = "px",width=3000, height=2000, res=350, pointsize=9)
par(mfrow=c(2,3))
scatterplots_metabo(module_color="yellow4",trait="Reticulocyte.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkslateblue",trait="Reticulocyte.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="lightyellow",trait="Reticulocyte.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="yellow")
dev.off()
png(file="TopModule_IntramodularStrength_Blood_Whitebloodcount.png", units = "px",width=3000, height=2000, res=350, pointsize=9)
par(mfrow=c(3,5))
scatterplots_metabo(module_color="bisque4",trait="Whitebloodcount.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="royalblue",trait="Whitebloodcount.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="lightpink4",trait="Whitebloodcount.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="lightcyan",trait="Whitebloodcount.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="blue")
scatterplots_metabo(module_color="mediumorchid",trait="Whitebloodcount.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="firebrick4",trait="Whitebloodcount.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkgreen",trait="Whitebloodcount.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkslateblue",trait="Whitebloodcount.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkviolet",trait="Whitebloodcount.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="skyblue1",trait="Whitebloodcount.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="darkorange2",trait="Whitebloodcount.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="lightcoral",trait="Whitebloodcount.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="antiquewhite4",trait="Whitebloodcount.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="skyblue",trait="Whitebloodcount.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
dev.off()
png(file="TopModule_IntramodularStrength_Blood_Bilirubin.png", units = "px",width=3000, height=2000, res=350, pointsize=9)
par(mfrow=c(2,2))
scatterplots_metabo(module_color="yellow4",trait="Bilirubin.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="lightcyan",trait="Bilirubin.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
scatterplots_metabo(module_color="skyblue1",trait="Bilirubin.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="")
dev.off()
png(file="TopModule_IntramodularStrength_Blood_LDH.png", units = "px",width=3000, height=2000, res=350, pointsize=9)
par(mfrow=c(1,2))
scatterplots_metabo(module_color="lightcyan",trait="LDH.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="blue")
scatterplots_metabo(module_color="steelblue",trait="LDH.beta",metabolite_significance=metabo_bloodtrait_cor_final, module_membership_cor=metabo_modulemembership_cor_final,alt_color="blue")
dev.off()